Fit trend models

Intro

This scripts fits models to gridded predictions of biomass trends as a function of climate trends (from climate-agnostic sdm, “04-fit-sdms-random.qmd”)

Load packages & source functions

# Load libraries, install if needed
pkgs <- c("tidyverse", "readxl", "tidylog", "RCurl", "devtools", "rMR", "ggridges",
          "kableExtra", "viridis", "RColorBrewer", "here", "sdmTMBextra") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Packages not on CRAN or dev version
# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)
library(sdmTMB)

# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")

# Set path
home <- here::here()

Read and plot data

# Read from GH?
d <- read_csv(paste0(home, "/data/clean/trends.csv")) |> 
  #mutate(quarter_f = as.factor(quarter)) |> 
  filter(quarter == 4) # probably there are ways to model this but for now we focus on Q4 because it's warmer and with less oxygen

Rows: 167822 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): id, species, life_stage, group
dbl (27): X, Y, quarter, temp_slope, biom_slope, oxy_slope, phi_slope, upr_t...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 87,036 rows (52%), 80,786 rows remaining

Fit sdmTMB models of biomass trends as a function of phi trends

#coefs <- list()
pred_trend1 <- list()
#pred_trend2 <- list()

for(i in unique(d$group)){
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15) # Need specific meshes because not for all hauls did we record all species
  
  print(
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
      labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
    )
  
    ggsave(paste0(home, "/figures/supp/mesh/10_trend_phi_", i, ".pdf"), width = 17, height = 17, units = "cm")

    print(i)
    
    m <- sdmTMB(biom_slope ~ phi_slope_sc * mean_phi_ct + mean_log_biomass_ct,
      mesh = mesh,
      data = dd,
      family = gaussian(), 
      #family = student(df = 4), 
      spatial = "on",
      spatiotemporal = "off"
      )
    
    sanity(m) 
  
    # Plot residuals
    # samps <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200)
    # mcmc_res <- residuals(m, type = "mle-mcmc", mcmc_samples = samps)
    # 
    # print(
    #   ggplot(dd, aes(sample = mcmc_res)) +
    #     stat_qq() +
    #     stat_qq_line() +
    #     labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
    #     theme(aspect.ratio = 1)
    #   )

    # Quick and dirty residuals
    res <- residuals(m)
    print(
      ggplot(dd, aes(sample = res)) +
        stat_qq() +
        stat_qq_line() +
        labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
        theme(aspect.ratio = 1)
      )
    
    ggsave(paste0(home, "/figures/supp/qq-trend/qq_trend_phi_", i, ".pdf"),
           width = 11, height = 11, units = "cm")
  
    # Save coefficients??
    #coefs[[i]] <- tidy(m, conf.int = TRUE) |> mutate(group = i)
    
    # Now predict the biotic trend for high and low temperature with a sd change in the climate trend
    low_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.05))
    high_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.95))
    mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
  
    # Predict
    nd <- data.frame(phi_slope_sc = c(-1, 0, -1, 0),
                     mean_phi_ct = c(low_mean_clim, low_mean_clim, high_mean_clim, high_mean_clim),
                     mean_log_biomass_ct = rep(mean_log_biomass_ct, 4))
    
    nsim <- 500
    pred <- predict(m, newdata = nd, nsim = nsim, re_form = NA)
    
    # one column is on iteration, one row is one nd row!
    pred <- data.frame(
      mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim), # see head(nd)
      est_change = c(pred[1, ], pred[3, ]), # see head(nd), these are declines in phi given low and high mean clim
      est_0 = c(pred[2, ], pred[4, ])) # see head(nd), these are the baseline average phi given low and high mean clim
    
    pred <- pred |> 
      mutate(delta_bio_trend = est_change - est_0, # if bio trend is positive, bio increases with declines in phi
             mean_clim = ifelse(mean_phi_ct == low_mean_clim, 
                                "low_phi", "high_phi"),
             group = i)
    
    pred_trend1[[i]] <- pred
    
    # Now make a prediction across a range of trends
    # Again I'm predicting for the range of the specific model
    # nd3 <- data.frame(expand.grid(phi_slope_sc = seq(quantile(dd$phi_slope_sc, probs = 0.05),
    #                                                quantile(dd$phi_slope_sc, probs = 0.95),
    #                                                length.out = 50),
    #                               mean_phi_ct = c(low_mean_clim, high_mean_clim))) |> 
    #   mutate(mean_log_biomass_ct = mean_log_biomass_ct)
    # 
    # pred3 <- predict(m, newdata = nd3, se_fit = TRUE, re_form = NA)
    # 
    # pred3 <- pred3 |> 
    #   mutate(lwr = est - 1.96*est_se,
    #          upr = est + 1.96*est_se, 
    #          group = i)
    # 
    # pred_trend2[[i]] <- pred3
    # 
    # ggplot(pred3, aes(phi_slope_sc, est, color = factor(mean_phi_ct))) + geom_line()
  
}

filter: removed 65,024 rows (80%), 15,762 rows remaining

[1] "Cod_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 2 rows containing non-finite values (`stat_qq()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,989 rows (89%), 8,797 rows remaining

[1] "Flounder_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 8 rows containing non-finite values (`stat_qq()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,339 rows (81%), 15,447 rows remaining

[1] "Flounder_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 31 rows containing non-finite values (`stat_qq()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,083 rows (81%), 15,703 rows remaining

[1] "Plaice_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 16 rows containing non-finite values (`stat_qq()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,159 rows (81%), 15,627 rows remaining

[1] "Plaice_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 18 rows containing non-finite values (`stat_qq()`).

Warning: Removed 18 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 18 rows containing non-finite values (`stat_qq()`).

Warning: Removed 18 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,336 rows (88%), 9,450 rows remaining

[1] "Cod_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

# phi_coefs <- dplyr::bind_rows(coefs) |> 
#   mutate(sig = "N",
#          sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig),
#          sig = ifelse(estimate > 0 & conf.low > 0, "Y", sig))

phi_pred_delta <- dplyr::bind_rows(pred_trend1) |> 
  separate(group, c("species", "life_stage"),
           sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))

mutate: no changes

# phi_pred <- dplyr::bind_rows(pred_trend2) |> 
#   separate(group, c("species", "life_stage"),
#            sep = "_", remove = FALSE) |> 
#   mutate(species = str_to_title(species),
#          life_stage = str_to_title(life_stage))

#write_csv(phi_coefs, paste0(home, "/output/phi_trend_coefs.csv"))
write_csv(phi_pred_delta, paste0(home, "/output/pred_phi_trend_delta.csv"))
#write_csv(phi_pred, "/output/pred_phi_trend.csv")

Fit sdmTMB models of biomass trends as a function of temperature trends, given mean temperature

pred_trend_temp <- list()
pred_trend_temp2 <- list()

for(i in unique(d$group)){
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15) # Need specific meshes because not for all hauls did we record all species
  
  print(
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
      labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
    )
  
    ggsave(paste0(home, "/figures/supp/mesh/10_trend_temp_", i, ".pdf"), width = 17, height = 17, units = "cm")

    print(i)
    
    m <- sdmTMB(biom_slope ~ temp_slope_sc * mean_temp_ct + mean_log_biomass_ct,
                mesh = mesh,
                data = dd,
                #family = gaussian(), 
                family = student(df = 7), 
                spatial = "on",
                spatiotemporal = "off")
    
    sanity(m)
    
    # Plot residuals
    # samps <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200)
    # mcmc_res <- residuals(m, type = "mle-mcmc", mcmc_samples = samps)
    # 
    # print(
    #   ggplot(dd, aes(sample = mcmc_res)) +
    #     stat_qq() +
    #     stat_qq_line() +
    #     labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
    #     theme(aspect.ratio = 1)
    #   )

    # Quick and dirty residuals
    res <- residuals(m)
    print(
      ggplot(dd, aes(sample = res)) +
        stat_qq() +
        stat_qq_line() +
        labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
        theme(aspect.ratio = 1)
      )
    
    ggsave(paste0(home, "/figures/supp/qq-trend/qq_trend_temp_", i, ".pdf"),
           width = 11, height = 11, units = "cm")
  
    # Predict
    low_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.05))
    high_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.95))
    mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
  
    nd <- data.frame(temp_slope_sc = c(0, 1, 0, 1),
                     mean_temp_ct = c(low_mean_clim, low_mean_clim, high_mean_clim, high_mean_clim),
                     mean_log_biomass_ct = rep(mean_log_biomass_ct, 4))
    
    nsim <- 500
    
    pred <- predict(m, newdata = nd, nsim = nsim, re_form = NA)
    
    pred <- data.frame(
      mean_temp_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
      est_0 = c(pred[1, ], pred[3, ]),
      est_change = c(pred[2, ], pred[4, ])) # change is increase
    
    pred <- pred |> 
      mutate(delta_bio_trend = est_change - est_0, # if delta_bio_trend is positive, biomass goes up with 1 increase in temp
             mean_clim = ifelse(mean_temp_ct == low_mean_clim, 
                                "low_temperature", "high_temperature"),
             group = i)
    
    pred_trend_temp[[i]] <- pred

}

filter: removed 65,024 rows (80%), 15,762 rows remaining

[1] "Cod_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,989 rows (89%), 8,797 rows remaining

[1] "Flounder_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,339 rows (81%), 15,447 rows remaining

[1] "Flounder_Juvenile"

Warning: The model may not have converged. Maximum final gradient:
0.102745737819191.

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✖ `b_j` gradient > 0.001

ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `b_j` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `b_j` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 13 rows containing non-finite values (`stat_qq()`).

Warning: Removed 13 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 13 rows containing non-finite values (`stat_qq()`).

Warning: Removed 13 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,083 rows (81%), 15,703 rows remaining

[1] "Plaice_Adult"

Warning: The model may not have converged. Maximum final gradient:
0.0100456365943273.

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✖ `b_j` gradient > 0.001

ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `b_j` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,159 rows (81%), 15,627 rows remaining

[1] "Plaice_Juvenile"

Warning: The model may not have converged. Maximum final gradient:
0.0150054043646151.

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✖ `b_j` gradient > 0.001

ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,336 rows (88%), 9,450 rows remaining

[1] "Cod_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

pred_trend_temp <- dplyr::bind_rows(pred_trend_temp) |> 
  separate(group, c("species", "life_stage"),
           sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))

mutate: no changes

write_csv(pred_trend_temp, paste0(home, "/output/pred_temp_trend_delta.csv"))

Fit sdmTMB models of biomass trends as a function of oxygen trends, given mean oxygen

pred_trend_oxy <- list()
pred_trend_oxy2 <- list()

for(i in unique(d$group)){
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15)
  
  print(
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
      labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
    )
  
    ggsave(paste0(home, "/figures/supp/mesh/10_trend_oxy_", i, ".pdf"), width = 17, height = 17, units = "cm")

    print(i)
    
    m <- sdmTMB(biom_slope ~ oxy_slope_sc * mean_oxy_ct + mean_log_biomass_ct,
                mesh = mesh,
                data = dd,
                family = gaussian(), 
                #family = student(df = 7),
                spatial = "on",
                spatiotemporal = "off")
  
    sanity(m) 
  
    # Plot residuals
    # samps <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200)
    # mcmc_res <- residuals(m, type = "mle-mcmc", mcmc_samples = samps)
    # 
    # print(
    #   ggplot(dd, aes(sample = mcmc_res)) +
    #     stat_qq() +
    #     stat_qq_line() +
    #     labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
    #     theme(aspect.ratio = 1)
    #   )
    
    # Quick and dirty residuals
    res <- residuals(m)
    print(
      ggplot(dd, aes(sample = res)) +
        stat_qq() +
        stat_qq_line() +
        labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
        theme(aspect.ratio = 1)
      )
    
    ggsave(paste0(home, "/figures/supp/qq-trend/qq_trend_oxy_", i, ".pdf"),
           width = 11, height = 11, units = "cm")
  
    # Predict
    low_mean_clim <- as.numeric(quantile(dd$mean_oxy_ct, prob = 0.05))
    high_mean_clim <- as.numeric(quantile(dd$mean_oxy_ct, prob = 0.95))
    mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
    
    nd <- data.frame(oxy_slope_sc = c(-1, 0, -1, 0),
                     mean_oxy_ct = c(low_mean_clim, low_mean_clim, high_mean_clim, high_mean_clim),
                     mean_log_biomass_ct = rep(mean_log_biomass_ct, 4))
    
    nsim <- 500
    
    pred <- predict(m, newdata = nd, nsim = nsim, re_form = NA)
    
    pred <- data.frame(
      mean_oxy_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
      est_change = c(pred[1, ], pred[3, ]),
      est_0 = c(pred[2, ], pred[4, ])
      )
    
    pred <- pred |> 
      mutate(delta_bio_trend = est_change - est_0,
             mean_clim = ifelse(mean_oxy_ct == low_mean_clim, 
                                "low_oxygen", "high_oxygen"),
             group = i)
    
    pred_trend_oxy[[i]] <- pred

}

filter: removed 65,024 rows (80%), 15,762 rows remaining

[1] "Cod_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 2 rows containing non-finite values (`stat_qq()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,989 rows (89%), 8,797 rows remaining

[1] "Flounder_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 8 rows containing non-finite values (`stat_qq()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,339 rows (81%), 15,447 rows remaining

[1] "Flounder_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 31 rows containing non-finite values (`stat_qq()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,083 rows (81%), 15,703 rows remaining

[1] "Plaice_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 16 rows containing non-finite values (`stat_qq()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,159 rows (81%), 15,627 rows remaining

[1] "Plaice_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 18 rows containing non-finite values (`stat_qq()`).

Warning: Removed 18 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 18 rows containing non-finite values (`stat_qq()`).

Warning: Removed 18 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,336 rows (88%), 9,450 rows remaining

[1] "Cod_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

pred_trend_oxy <- dplyr::bind_rows(pred_trend_oxy) |> 
  separate(group, c("species", "life_stage"),
           sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))

mutate: no changes

write_csv(pred_trend_oxy, paste0(home, "/output/pred_oxy_trend_delta.csv"))

What is the effect of temperature given oxygen? first check correlation

d |> group_by(group) |> summarise(cor = cor(temp_slope_sc, oxy_slope_sc))

group_by: one grouping variable (group)

summarise: now 6 rows and 2 columns, ungrouped

# A tibble: 6 × 2
  group               cor
  <chr>             <dbl>
1 Cod_Adult         0.240
2 Cod_Juvenile      0.182
3 Flounder_Adult    0.399
4 Flounder_Juvenile 0.185
5 Plaice_Adult      0.177
6 Plaice_Juvenile   0.165

Fit sdmTMB models of biomass trends as a function of temperature trends, given mean oxygen

pred_trend_temp_o <- list()
pred_trend_temp2_o <- list()

for(i in unique(d$group)){
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15) # Need specific meshes because not for all hauls did we record all species
  
  print(
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
      labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
    )
  
    ggsave(paste0(home, "/figures/supp/mesh/10_trend_temp_o_", i, ".pdf"), width = 17, height = 17, units = "cm")

    print(i)
    
    m <- sdmTMB(biom_slope ~ temp_slope_sc * mean_oxy_ct + mean_log_biomass_ct,
                mesh = mesh,
                data = dd,
                family = gaussian(), 
                #family = student(df = 7), 
                spatial = "on",
                spatiotemporal = "off")
    
    sanity(m)
    
    # Plot residuals
    # samps <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200)
    # mcmc_res <- residuals(m, type = "mle-mcmc", mcmc_samples = samps)
    # 
    # print(
    #   ggplot(dd, aes(sample = mcmc_res)) +
    #     stat_qq() +
    #     stat_qq_line() +
    #     labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
    #     theme(aspect.ratio = 1)
    #   )
  
    # Quick and dirty residuals
    res <- residuals(m)
    print(
      ggplot(dd, aes(sample = res)) +
        stat_qq() +
        stat_qq_line() +
        labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
        theme(aspect.ratio = 1)
      )
    
    ggsave(paste0(home, "/figures/supp/qq-trend/qq_trend_temp_o_", i, ".pdf"),
           width = 11, height = 11, units = "cm")
  
    # Predict
    low_mean_clim <- as.numeric(quantile(dd$mean_oxy_ct, prob = 0.05))
    high_mean_clim <- as.numeric(quantile(dd$mean_oxy_ct, prob = 0.95))
    mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
  
    nd <- data.frame(temp_slope_sc = c(0, 1, 0, 1),
                     mean_oxy_ct = c(low_mean_clim, low_mean_clim, high_mean_clim, high_mean_clim),
                     mean_log_biomass_ct = rep(mean_log_biomass_ct, 4))
    
    nsim <- 500
    
    pred <- predict(m, newdata = nd, nsim = nsim, re_form = NA)
    
    pred <- data.frame(
      mean_oxy_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
      est_0 = c(pred[1, ], pred[3, ]),
      est_change = c(pred[2, ], pred[4, ])) # change is increase in temp
    
    pred <- pred |> 
      mutate(delta_bio_trend = est_change - est_0, # if delta_bio_vel is positive, biomass goes up with 1 increase in temp
             mean_clim = ifelse(mean_oxy_ct == low_mean_clim, 
                                "low_oxygen", "high_oxygen"),
             group = i)
    
    pred_trend_temp_o[[i]] <- pred

}

filter: removed 65,024 rows (80%), 15,762 rows remaining

[1] "Cod_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 2 rows containing non-finite values (`stat_qq()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq()`).

Warning: Removed 2 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,989 rows (89%), 8,797 rows remaining

[1] "Flounder_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 8 rows containing non-finite values (`stat_qq()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq()`).

Warning: Removed 8 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,339 rows (81%), 15,447 rows remaining

[1] "Flounder_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 31 rows containing non-finite values (`stat_qq()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq()`).

Warning: Removed 31 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,083 rows (81%), 15,703 rows remaining

[1] "Plaice_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 16 rows containing non-finite values (`stat_qq()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq()`).

Warning: Removed 16 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 65,159 rows (81%), 15,627 rows remaining

[1] "Plaice_Juvenile"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Warning: Removed 19 rows containing non-finite values (`stat_qq()`).

Warning: Removed 19 rows containing non-finite values (`stat_qq_line()`).

Warning: Removed 19 rows containing non-finite values (`stat_qq()`).

Warning: Removed 19 rows containing non-finite values (`stat_qq_line()`).

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

filter: removed 71,336 rows (88%), 9,450 rows remaining

[1] "Cod_Adult"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

mutate: new variable 'delta_bio_trend' (double) with 1,000 unique values and 0% NA

        new variable 'mean_clim' (character) with 2 unique values and 0% NA

        new variable 'group' (character) with one unique value and 0% NA

pred_trend_temp_o <- dplyr::bind_rows(pred_trend_temp_o) |> 
  separate(group, c("species", "life_stage"),
           sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))

mutate: no changes

write_csv(pred_trend_temp, paste0(home, "/output/pred_temp_o_trend_delta.csv"))

knitr::knit_exit()